﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using Microsoft.Xna.Framework;

namespace COC.Util
{
	public class SunCalculator
	{
		private static float m_dbSunElevation = MathHelper.Pi / 4;
		private static float m_dbSunDistance = 150000000000;
		private static float m_dbSunHeading = -MathHelper.Pi / 4;

		public static float SunElevation
		{
			get { return m_dbSunElevation; }
			set { m_dbSunElevation = value; }
		}

		public static float SunDistance
		{
			get { return m_dbSunDistance; }
			set { m_dbSunDistance = value; }
		}

		public static float SunHeading
		{
			get { return m_dbSunHeading; }
			set { m_dbSunHeading = value; }
		}

		public static Vector3 GetGeocentricPosition(
			DateTime utcDateTime,
			float dbCameraLat, float dbCameraLon)
		{
			// Sun synched with time and date
			float JD = getJulianDay(utcDateTime);
			float T = (JD - 2451545.0f) / 36525f; // number of Julian centuries since Jan 1, 2000, 12 UT

			float k = MathHelper.Pi / 180.0f;
			float M = 357.52910f + 35999.05030f * T - 0.0001559f * T * T - 0.00000048f * T * T * T; // mean anomaly, degree
			float L0 = 280.46645f + 36000.76983f * T + 0.0003032f * T * T; // mean longitude, degree
			float DL = (float)((1.914600f - 0.004817f * T - 0.000014f * T * T) * Math.Sin(k * M) + (0.019993f - 0.000101f * T) * Math.Sin(k * 2 * M) + 0.000290f * Math.Sin(k * 3 * M));
			float L = L0 + DL; // true longitude, degree
			L = L % 360;
			//		obliquity eps of ecliptic:
			float eps = 23.0f + 26.0f / 60.0f + 21.448f / 3600.0f - (46.8150f * T + 0.00059f * T * T - 0.001813f * T * T * T) / 3600;

			float X = (float)Math.Cos(k * L);
			float Y = (float)(Math.Cos(k * eps) * Math.Sin(k * L));
			float Z = (float)(Math.Sin(k * eps) * Math.Sin(k * L));

			float R = (float)Math.Sqrt(1.0 - Z * Z);

			float dec = (float)((180 / MathHelper.Pi) * Math.Atan(Z / R)); // in degrees
			float RA = (float)((24 / MathHelper.Pi) * Math.Atan(Y / (X + R))); // in hours

			float theta0 = 280.46061837f + 360.98564736629f * (JD - 2451545.0f) + 0.000387933f * T * T - T * T * T / 38710000.0f; // Greenwich Sidereal Time

			theta0 = theta0 % 360;
			RA *= 15; // convert from hours to degrees
			float tau = theta0 - RA;

			float latRadians = MathHelper.ToRadians(-dec);
			float lonRadians = MathHelper.ToRadians(-(tau - 180));
			float radCosLat = (float)Math.Cos(latRadians);

			Vector3 v3Pos = new Vector3(
				radCosLat * (float)Math.Cos(lonRadians),
				radCosLat * (float)Math.Sin(lonRadians),
				(float)Math.Sin(latRadians));


			return v3Pos;
		}

		private static float getJulianDay(System.DateTime gregorianDate)
		{
			float d = gregorianDate.Day +
				(float)gregorianDate.Hour / 24.0f +
				(float)gregorianDate.Minute / (24.0f * 60.0f) +
				(float)gregorianDate.Second / (24.0f * 60.0f * 60.0f) +
				(float)gregorianDate.Millisecond / (24.0f * 60.0f * 60.0f * 1000.0f);

			int m = (gregorianDate.Month < 3 ? gregorianDate.Month + 12 : gregorianDate.Month);
			int j = (gregorianDate.Month < 3 ? gregorianDate.Year - 1 : gregorianDate.Year);

			int c = 2 - (int)Math.Floor((float)(j / 100)) + (int)Math.Floor((float)(j / 400));
			return (int)Math.Floor((float)(1461 * (j + 4716)) / 4) + (int)Math.Floor((float)(153 * (m + 1) / 5)) + d + c - 1524.5f;
		}
	}
}
